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The shear stress autocorrelation function has been studied recently by 

molecular dynamics simulation using the q~ n potential for very large n. The 

results are analyzed and interpreted here by comparing them to the shear 

stress response function for hard spheres. It is shown that the hard sphere 

response function has a singular contribution and that this is reproduced 

accurately by the simulations for large n. A simple model for the stress 

autocorrelation function at finite n is proposed, based on the required hard 

sphere limiting form. 



I. INTRODUCTION 

The qualitative features of simple atomic fluids are often modelled using the hard sphere 
pair potential: zero for q > a and infinite for q < o. It is well-known that the thermody- 
namic and structural features of the hard sphere fluid are represented continuously by the 
soft sphere q~ n potential in the limit n — > oo. In contrast, the dynamical properties for 
hard and soft sphere fluids are different since the collision time (time required for finite 
momentum transfer) is strictly zero in the former case and of order 1/n in the latter case. 
As shown below, this has a significant effect on the short time behavior of time correla- 
tion functions. Recently, Powles and co-workers |l]-[3] have investigated the time correlation 
functions determining the shear and bulk viscosity for the q~ n fluid at large values of n (up 
to n = 288). The initial values of these functions diverges oc n, and the simulation results 
for the normalized correlation functions appear to scale with the single characteristic time 
r n oc 1/n. These results are puzzling when considered as an approach to the hard sphere 
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limit. The initial values, i.e. the equal time stress tensor fluctuations, are related to elastic 
constants which should be finite for the hard sphere fluid. Also, the only time scale in the 
problem r n vanishes for large n leaving no finite time dynamics for the hard spheres fluid. In 
this picture, the finite Green-Kubo shear viscosity results from the product of the diverging 
initial condition and a vanishing correlation time. Clearly, the relationship of properties at 
large n to those at n = oo must be more complex. The objective here is to resolve these 
paradoxes by calculating the corresponding properties for the hard sphere fluid directly to 
allow comparison with those at large n. Attention is restricted to the response function for 
shear viscosity, but similar results apply for the bulk viscosity response function as well. 

Standard methods of linear response lead to Green-Kubo expressions for the shear viscos- 
ity as the time integral of a shear response function. For soft spheres this response function 
is proportional to the autocorrelation function for the microscopic stress tensor. The form of 
the stress tensor is identified from the micrscopic local conservation law for the momentum 
density. Since it explicitly involves the interparticle force, a divergence for large n can be 
anticipated. Application of the same linear response methods for the hard sphere fluid lead 
to several differences: 

• The response function is not simply the stress tensor autocorrelation function, as for 
soft spheres, but has an additional singular contribution proportional to S(t). 

• The coefficient of the singular contribution can be calculated exactly and is simply 
related to the limiting form of the stress fluctuations for the soft sphere fluid. 

• The form of the microscopic stress tensor for hard spheres is quite different from that 
for soft spheres, and the initial value of the stress tensor autocorrelation function is 
finite. 

• The hard sphere stress tensor autocorrelation function has a characteristic time scale 
of the mean free time, which is the same as that for the soft sphere fluid at large n. 
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The singular part of the hard sphere correlation function evolves continuously from the soft 
sphere correlation function as n —>■ oo, and is closely related to the divergent initial value. 
The results of Powles et al. are shown to be an accurate measurement of this singular 
part. However, the resolution of this component in the simulation tends to mask the finite 
remainder that persists on the mean free path time scale, whose existence is clear from 
the hard sphere calculations and which is essential for the correct density dependence of 
transport properties. An alternative representation of the simulation results is proposed to 
more closely demonstrate the approach to the hard sphere limit. 

The relevant features of the dynamics and statistical mechanics of hard spheres are in- 
troduced in the next section to provide the necessary definitions and to contrast with the 
corresponding features for soft spheres. In Section III the results for linear response to a 
transverse velocity field perturbation are quoted, and the Green-Kubo expression is identi- 
fied. The singular contribution is evaluated exactly and given a representation appropriate 
for comparison with the MD simulation results. The remainder, determined from the stress 
tensor autocorrelation function for hard spheres is evaluated in an Enskog approximation. 
The limiting form for soft spheres approaching the hard sphere limit is discussed on the 
basis of these results. Some final comments are offered in Section IV. 

It is a pleasure to dedicate this brief contribution to Keith Gubbins - friend, mentor, and 
colleague - who has played such an important role in unfolding the mysteries of transport 
in simple and complex fluids. 

II. STATISTICAL MECHANICS OF HARD SPHERES 

Any interatomic potential with discontinuities (e.g., hard sphere, square well) implies an 
infinite force at the point of discontinuity. Consequently, the usual formulation of Newtonian 
dynamics in terms of forces alone is no longer appropriate. In this section the dynamics and 
statistical mechanics for a system of iV smooth hard spheres with diameter a is reviewed 
briefly §. 
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The dynamics consists of free streaming until a given pair i, j is in contact, at which point 
the relative velocity of that pair changes instantaneously according to the elastic collision 
rule 



g 



gij - 2S : (g ii • a) 



(1) 



Here = Vj — Vj is the relative velocity for particles i and j. The state of the sys- 
tem at time t is completely characterized by the positions and velocities of all spheres at 
that time and is represented by a point in the associated 6N dimensional phase space, 
T t = {qi(t), .., qjv(t), Vi(i), .., Vjv(t)}. The sequence of free streaming and binary collisions 
determines uniquely these positions and velocities of the spheres at time t for given initial 
conditions. A more complete notation expressing this dependence on initial conditions is 
r t (IV). Observables of interest are represented by the phase functions A(T t (T)), and their 
average for given statistical initial data at t = is defined by 



where p(T) is the probability density or ensemble for the initial state, normalized to unity. 
An equivalent representation of this average is obtained by changing variables to integrate 
over Tt rather than over T expressed as a function of T t , i.e. T = iy 1 {Tt)- This change 
of variables in (§) allows the time dependence to be expressed in terms of the probability 
density 



where the probability density at time t is defined by p(T,t) = p{T~[ (T)). The Jacobian of 
the transformation is unity. All of the above is the same for both hard and soft spheres. 
The equilibrium time correlation functions for observables A and B are defined by 




(2) 




(3) 



C AB (t) =< 5A(t)8B > 




dT Pe (T)SA(t)SB, 6X 



= X— < X > e . 



(4) 



where A (t) = A(T t (To)). The Gibbs distribution given by 
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p.(T) = W({fl«})n^(«0, = (^) 3/2e_/3m,;2/2 ( 5 ) 



is a stationary state for the hard sphere dynamics. The "overlap" function W({qij}) is 
defined to vanish if any pair overlap, < a, and is normalized to unity when integrated 
over all configuration space. The equilibrium correlation function is stationary 

< 8A(t)6B > e =< 5A8B(-t) > e (6) 

where B(—t) = B{T-t(T)) has the time reversed dynamics. 

For practical purposes it is useful to identify the generators L± and L for the two 
representations, defined for t > by 

A(±t) = e ±L ^A(T), p(r, t) = e- Tt p(T). (7) 

These are not the usual generators of Hamilton's equations for continuous forces and are 
somewhat more complex due to the singular nature of hard spheres. They are given by 

N N N N N N 

i=l i=l jj^i i=l 1=1 j^i 

The first terms on the right sides generate free streaming while the second terms describe 
velocity changes. The three binary collision operators, T±(i,j) and T(i,j), for particles i 
and j are given by 

T±(iJ) = W J dtt e( Tgii -5 : )(g i ,-5 ; )5(q ij - a)^ - 1), (9) 

T(i,j) =a 2 J dVl e(g ii -ff)(gi i -ff)[5(q y - <r)% - %y + <x)], (10) 

where dfl denotes the solid angle integration for the unit vector &, r is the relative 
position vector of the two particles, and the operator by is a substitution operator, 
bijF(gij)=F(bijgij), which changes the relative velocity gij = v, — Vj into its scattered 
velocity according to (|3]) fr^gy = gTJ. For continuous potentials all generators are the same 
L± = L but for hard spheres all three generators clearly are different. 



The local microscopic conservation laws follow from the above dynamics by differentiating 
the local conserved densities of mass, energy, and momentum with respect to time. The 
relevant conservation law for the discussion here is that for momentum density 

d tPa (r, ±t) + djT ±a p(r, ±t) = 0, (11) 

N 

P( r ) = Yl Pi< ^ r ~ d j T ±*p( r ) = TL±p a {r) (12) 

i=l 

The subscript ± distinguishes the forms of the momentum flux for the case t > and t < 0. 
The detailed expression for r±ij(r,t) is not required below, only its volume integral 

T± a p = J drT± af3 (r) 

N 1 N 

= Y m Vi, a Vj,p + -cr 3 ^2 / dQ Q(Tgiyd)m(g ir d) 2 a a ap5(q ij - a) (13) 
i=i i^j J 

This form for the hard sphere momentum flux(es) can be contrasted with that for the soft 

sphere potential v(q) = e (cr/q) n 

Tap = ^2 mv i>a Vifi + - Y 9ija Qijpne I — J (14) 

which can be written for large n as 

N l N 

T a p = ^2 mVi^Vi^ + n- / dQ a a dp5(qij - a) (15) 
i=\ i^j J 

While flTB] ) and (|IJJ) are similar, both involving a surface integral for pairs of particles at 

contact, there are two significant differences. The hard sphere stress tensor is momentum 

dependent and the soft sphere stress tensor is proportional to n, diverging in the large n 

limit. 



III. LINEAR RESPONSE 



The results of the last section now allow a discussion of the response functions for the hard 
sphere fluid. Attention here is restricted to the response to an initial transverse perturbation 
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of the velocity field, u y (x,0) = Asin(x/X). The velocity field for A much larger than the 
mean free path and for times long compared to the mean free time then obeys the diffusion 
equation 

d t u y (x, t) + (ji/n e m) V 2 u y (x, t) = (16) 

where r\ is the shear viscosity and n e is the equilibrium density . The derivation of this 
result is straightforward using the standard methods of linear response and the Liouville 
dynamics defined in the last section. The details will be given elsewhere || and only the 
results discussed here. The shear viscosity is given by the Green-Kubo expression 



V = P 

where 



1 f°° 1 

y (T +X yA xy ) e + dt— (T +xy (t)T_ xy ) e 



(17) 



A xy = / drxp y (r) = ^2q ix mv iy (18) 

For soft spheres T +xy = T- xy = % y and (T xy A xy ) e = since T xy is independent of the particle 
velocities while A xy is an odd function of the velocities. Thus, the usual Green-Kubo form 
for the shear viscosity in terms of the time integral of the stress tensor autocorrelation 
function is regained 



v ss = P 



poo 1 

J o dt-{T xy (t)T xy ) e (19) 



In contrast, (T +xy A xy ) e does not vanish for hard spheres. The equivalence of fll7D and (|19D 
in the limit of n — > oo requires 

lim i (% y (t) % v ) e = ~ (T +xy A xy ) e 5(t) + i (T +xy (t) T„ xy ) e (20) 
n — >oo y y y 

This is the primary observation of this work regarding the relationship of shear response for 
soft and hard spheres. The hard and soft sphere stress autocorrelation functions differ by a 
singular term. 

To explore this relationship, two results are noted. First, for hard spheres the coefficient 
of the singular term can be calculated exactly to get || 



- - ( WW, = -a ^— j (p - p k ) . 



(21) 



Also, for soft spheres the exact initial condition is 
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y (Txy (0) %cy) e = ^^oo — -jj 

Here (3 = 1/ksT, p k = n e //3, p is the pressure 
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P k + - (n - 3) (p - p*) 



(22) 



P = P I 1 + g^eX 



(23) 



and x is the equilibrium pair correlation function evaluated at = a. Equation (p2|) shows 
the divergence of the high frequency shear modulus, Goo, anti — ► oo. Furthermore, the 
dominant short time behavior of (T xy (t) T xy ) for soft spheres with large n is 
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\%cy{t)%cy) e (T xy (fS) T xy ) ( 

with the characteristic short time scale 



+ order i 4 



(24) 



ay/ (3m T\ 



n 



n 



(25) 



Comparison of (pT|) and (|22|) shows the relationship 

4r n 1 



lim 

n— >oo ^/7T V/ 



(^zj/ (0) Txy) £ ( y (T+xy^-xy) e ^ 



(26) 



Note that the left side is finite in this limit due to the factor of r n ~ 1/n. 

These results provide an appropriate means to discuss and interpret the simulation results 
of reference 0. The stress tensor autocorrelation function is written as 



(%y (t) Txy) e — (0) 7iy) e C n (t) 



(27) 



and the simulation data is reported for the normalized correlation function C n (t) (a subscript 
n has been included since the simulations are performed for n = 18, 36, 72, 144, and 288). It 
is found that the results are globally well-fit by 

t 



[Cn(t)i 



sech a 



(28) 



where a = 7r 3 / 2 /4 « 1.4 (in reference || a = v^2)- This is the puzzling result mentioned in 
the Introduction above. The coefficient of C n (t) in (|27| ) diverges for large n while [C(t)] sim is 
a "universal" function of time with a vanishing time scale. To understand this paradoxical 
behavior combine (|26| ) - (^) to get, for large n 



V < T +- A ^> e ) if sech (4) 

(-y{T +xy A xy )^5 n {t) (29) 



In the second line 5 n (t) is a representation of the delta function as n — > oo defined by 

5 n (t) = — sech (a—) (30) 



vrr„ \ t. 



It is now clear that the MD simulation has measured the singular contribution on the 
right side of (|20|). This reinterpretation of the MD results shows a remarkable confirmation 
of the hard sphere singular contribution evolving continuously from the soft sphere stress 
autocorrelation function. It also demonstrates the role of the diverging shear modulus in 
constructing the amplitude of 5 n (t). Finally, the paradox of the vanishing time scale is now 
understood as well. The second term on the right side of ( |20"D is expected to have a time 
scale of the order of the mean free time (see below), but has not been well- resolved by the 
simulations. The normalization to define C(t), appropriate to identify this singular term, 
makes observation of the second term difficult since it then vanishes as n" 1 . Still, such 
deviations from the fit (p8|) are observed in the simulations, as noted below. 

The second term in (p0[) is the hard sphere stress autocorrelation function. It can be 
evaluated approximately from Enskog kinetic theory. The details will be given elsewhere, 
and only the result quoted here 

\ (T +xy (t) T_ xy ) e - /TV (l + l^f) 2 e-* /T »* (31) 
where T m f p is the mean free time 

5 p k 



imfp 
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r ^n. (32) 



The pressures for the hard and soft sphere fluids are continuously related for large n so 
both the non-singular part in fl3~T|) and the mean free path should be similar in both cases. 
It is easily verified that use of ( |31~D and ( [21] ) in the Green-Kubo expression ( |TTD yields the 
exact Enskog shear viscosity. Both the singular term and the stress auto correlation function 
contribute at intermediate velocities. In the low density limit the singular contribution to 
the viscosity vanishes and only the stress autocorrelation function contributes. Also, note 
that the value at t = for the hard sphere stress autocorrelation function is finite, in contrast 
to that for the soft sphere correlation function. The implications of this for elastic constants 
will be discussed elsewhere. 

The hard sphere results suggest how to write an approximate stress autocorrelation 
function for the soft sphere potential at large n. There should be two contributions with 
two different time scales. The first is that which has been determined from the simulations 
and which yields the singular contribution in the hard sphere limit. It also gives entirely 
the exact short time behavior. The second term should vanish at short times, but tend 
toward the form ( |3~T| ) for t » r n . A simple approximation consistent with the above exact 
properties and which has the proper hard sphere limit is 

I<r,(.)T,).*I Go *e*( a l) ( 1 + ?^)V'/-, (33) 

Equivalently, the normalized correlation function is 

C n it) « sech (®-^J + A n (t) (34) 

where A n (t) is the contribution from the non-singular part 

An (t) = X{-)£- ( 1 + |£^V e-t/w, (35) 



X V 5 Vk 

The function X(t) controls the crossover from the time scale r n to the mean free time scale. 
For consistency it must be even in time, vanish up through order t 2 , and approach unity for 
large times. One simple form with these properties is 

*M=(l^) 2 (36) 
10 
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FIG. 1. Plot of normalized stress autocorrelation function C(t) vs log t at packing fraction 0.45 



for n = 18, 76, and 288 
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FIG. 2. Same as Figl for A(t), non-singular contribution to C(t) 



For illustration, Figure 1 shows C n (t) as a function of t/r n for n = 18, 76, and 288 at the 
packing fraction irn e a 3 /6 = 0.45. This is quite similar to the simulation results in Figure 5 
of reference 0. The small variations with n for t/r n > 1 are due to A n (£). The latter is 
shown in Figure 2. 
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A more rational choice for the crossover function X(t) is obtained by relating it to the 
low density autocorrelation function 



X { — ) = e ' m fp 



C® (t) - sech (a-^j 



(37) 



The superscript (0) denotes the low density limit. For example, Tn^ is the time scale obtained 
from the expansion of the low density autocorrelation function to second order in t. With 
this choice the normalized stress autocorrelation function becomes 



C n (i) « sech ( a— ) + 



C {0) (t) - sech 



t 



r (o) 

In 



where the new time scale of the second term is 



24 



T T m f p T 



(0) 



mfp 



P 



p — p 
Pk 



5 Pk 



(38) 



(39) 



The calculation of (t) for all times is still a difficult problem, but one that has been 
reduced to two particle dynamics. The remaining terms are all expressed in terms of the 
pressure. It should be interesting to see how well ( pB]) predicts the density dependence of 
the viscosity for the soft sphere fluid. 



IV. DISCUSSION 

On purely physical grounds it is clear that the physical properties of the soft sphere fluid 
for large n should approach those of the hard sphere fluid, in general. This requires that 
the continuous Newtonian dynamics of the soft sphere fluid must generate the singularities 
that are inherent in the dynamics of the hard sphere fluid. The molecular dynamics simu- 
lation results of references [|l]-|3|] provide the first direct demonstration of how such singular 
properties can be extracted from competition between diverging fluctuations (the high fre- 
quency shear modulus) and a vanishing collision time. More generally, the autocorrelation 
functions for the soft sphere fluid have contributions on two time scales, r n and r m f p . For 
asymptotically short times only the former contributes. There is then a crossover to longer 
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times when the latter dominates. As the density increases, r m / p decreases and if n is not 
large the two times can become comparable and the dynamics is a composition of both. This 
is the case in Figure 2 for 1 ^ t/r n J$ 10. A reevaluation of the simulation data to extract 
the non-singular dynamics on the time scale r m f p is of some interest and could be compared 
directly to hard sphere simulations at t > 0. It appears that a recently proposed approxi- 
mation [0] to interpolate correlation functions between their exact short time dynamics and 
long time exponential decay does not work in this case due to the divergent initial value. 

Similar questions regarding the relationship of hard and soft sphere systems can be raised 
for the solid state. The elastic constants are usually related to the stress tensor fluctuations 
which diverge for the soft sphere fluid but are finite for the hard sphere fluid. To clarify 
this, the exact relations (p0[ ) and ( pB"D can be combined to give 

= v -(T +xy (t)T_ xy ) e (40) 

The right side is finite for t = + ([J^), which suggests that the proper elastic properties 
of the soft sphere fluid require subtraction of the short time dynamics with the diverging 
amplitude. This can be viewed as an annealing process on the time scale r n . 

The analysis of simulation data for thermal conductivity and bulk viscosity for the soft 
sphere fluid leads to similar results. The description will be given in a more detailed report 
on this topic elsewhere ||. 
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